!************************************************************************************************************
! This module contains the inner workings of the halo model, including
! substructure of various incarnations: (power law, delta function, and
! Gaussian sub-halo mass functions).
!
! Declare a tHalo type:
!
! type(tHalo) :: halo
!
! and set it up with new:
!
! call new(halo, [lps], [fracSub], [mu], [cDistribution])
!
! You can specify your own lps type, or it will use the default (see base/linearPower.F90).
! [fracSub] to specify the fraction of a halo's mass in substructure (default 0%)
! [mu] is the exponent in substructure mass function (default 0.9)
! [cDistribution] to .true. if you want the log-normal distribution in halo concentration parameters
!
! Now you can calculate the power spectrum using the halo model using calls like:
!
! call powerSpectrum(halo, k, a, answer)
!
! Also, you can calculate the contributions separately.
!
!************************************************************************************************************
module mHaloModel

        use mPrecision
        use mObject
        use mNu

        implicit none

        interface new
                module procedure new_halo
        end interface new

        interface write
                module procedure write_halo
        end interface write

        interface getUniverse
                module procedure getUniverse_halo
        end interface getUniverse

        interface getAvgDensity
                module procedure getAvgDensity_halo
        end interface getAvgDensity

        interface getH
                module procedure getH_halo
        end interface getH

        interface getOmega_m0
                module procedure getOmega_m0_halo
        end interface getOmega_m0

        interface nu_m
                module procedure nu_m_halo
        end interface nu_m

        interface m_nu
                module procedure m_nu_halo
        end interface m_nu

        interface getLPS
                module procedure getLPS_halo
        end interface getLPS

        interface getSpectralIndex
                module procedure getSpectralIndex_halo
        end interface getSpectralIndex

        interface getMass_nu
                module procedure getMass_nu_halo
        end interface getMass_nu

        interface getRun
                module procedure getRun_halo
        end interface getRun

        interface dnu
                module procedure dnu_halo
        end interface dnu

        real(D), parameter :: MLL = 1e3_D, &
                              MUL = 1e16_D, &
                              SUBHALO_UL_FRACTION = 0.01_D, &
                              ! lightest subhalo that can have a galaxy:
                              ! note: that I set M_UL_FOR_G to subUL:
!                             M_LL_FOR_G = 1e13_D &
!                             M_LL_FOR_G = 1e12_D &
!                             M_LL_FOR_G = 5e11_D &
!                             M_LL_FOR_G = 1e11_D &
!                             M_LL_FOR_G = 8e10_D &
!                             M_LL_FOR_G = 5e10_D &
!                             M_LL_FOR_G = 1e10_D &
!                             M_LL_FOR_G = 5e9_D &
!                             M_LL_FOR_G = 1e9_D &
                             M_LL_FOR_G = 1e8_D, &
!                             M_UL_FOR_G = 1e14_D, &
                             M_UL_FOR_G = 0.9e16_D, &
                             SIGMA_LN_C = 0.3_D, &         ! for c
                             SIGMA_LN_C2 = SIGMA_LN_C ** 2 ! distribution
#ifdef M_VIR
        ! New fits from White (2002):
        real(D), parameter :: SHETH_TORMEN_A = 0.79_D, &
                              SHETH_TORMEN_P = 0.32_D, &
            ! Norm for p = 0.32:
                SHETH_TORMEN_NORM = 0.1202086379800250895358667946943854_D
        ! Old parameters:
!        real(D), parameter :: SHETH_TORMEN_A = 0.707_D, &
!        real(D), parameter :: SHETH_TORMEN_A = 0.75_D, &
!                              SHETH_TORMEN_P = 0.3_D, &
            ! Norm for p = 0.3:
!                SHETH_TORMEN_NORM = 0.1285326740593642818863074431355724_D &
#else ! not M_VIR
        ! New fits from White (2002):
        real(D), parameter :: SHETH_TORMEN_A = 0.67_D, &
                              SHETH_TORMEN_P = 0.33_D, &
                SHETH_TORMEN_NORM = 0.1157617500173089837953450055285082_D, &
                              ! Mean overdensity of a halo relative to
                              ! background density
                              DELTA_H = 180.0_D, &
                DELTA_H_TO_ONE_THIRD = 5.646216173286170946465316362087775_D
#endif ! not M_VIR

#ifdef DEV
!        real(D), parameter :: EPS = 1e-6_D
!        real(D), parameter :: EPS = 1e-5_D
!        real(D), parameter :: EPS = 1e-4_D
        real(D), parameter :: EPS = 1e-3_D
!        real(D), parameter :: EPS = 1e-2_D
#else ! DEV
!        real(D), parameter :: EPS = real( epsilon(1.0_R4), D )
!        real(D), parameter :: EPS = 1e-6_D
        real(D), parameter :: EPS = 1e-5_D
!        real(D), parameter :: EPS = 1e-4_D
!        real(D), parameter :: EPS = 1e-3_D
!        real(D), parameter :: EPS = 1e-2_D
!        real(D), parameter :: EPS = 1e-1_D
#endif ! DEV

        private
        public :: new, write
        public :: getUniverse, getH, getAvgDensity
        public :: getFracSub, setFracSub
        public :: getM_star
        public :: DELTA_H
        public :: normalizedForA, setNormalizedForA
        public :: cDistribution, setCDistribution
        public :: normalizeBias, biasNorm, setBiasNorm
        public :: getLPS
        public :: getSpectralIndex, getMass_nu, getRun
        public :: subLL, subUL
        public :: SIGMA_LN_C, SIGMA_LN_C2, EPS
        public :: m_nu, nu_m, nuUL, nuLL
        public :: refitNu
        public :: MLL, MUL, SUBHALO_UL_FRACTION
        public :: haloProfile, NFW, cParam, dimensionless_mProfile, STf_Dnu
        public :: N_0, subMassFunctionHm
        public :: m_r_s, m_R_vir
        public :: PSf_dnu, PS_b1, ST_b1
        public :: R_M_180, m_r_180
        public :: N_gal, N_sub, M_LL_FOR_G, M_UL_FOR_G, alpha
        public :: checkBiasNorm
        public :: dnu
        public :: powerSpectrum, P_1h, P_2h, P_1c, P_sc, P_2c
!***************************************************************************
contains

        subroutine new_halo(halo, lps, fracSub, mu, cDistribution)

                use mLPS

                type(tHalo), intent(out) :: halo
                type(tLPS), intent(in), optional :: lps
                real(D), intent(in), optional :: fracSub,mu
                logical, intent(in), optional :: cDistribution

                real(D), parameter :: DEFAULT_FRAC_SUB = 0.0_D, &
                                      DEFAULT_MU = 0.9_D

                if( present(lps) ) then
                        call new(halo%nuFit, lps, a = 1.0_D)
                else
                        call new(halo%nuFit, a = 1.0_D)
                end if

                if( present(fracSub) ) then
                        halo%fracSub = fracSub
                else
                        halo%fracSub = DEFAULT_FRAC_SUB
                end if

                if( present(mu) ) then
                        halo%mu = mu
                else
                        halo%mu = DEFAULT_MU
                end if

                if( present(cDistribution) ) then
                        halo%cDistribution = cDistribution
                else
                        halo%cDistribution = .false.
                end if

!DEREK                halo%biasNorm = 1.0_D
!                halo%normalizedForA = -1.0_D
                call normalizeBias(halo, 1.0_D)

                halo%m_star = nu_m(halo, 1.0_D, 1.0_D)
                call setupSubhaloLimits(halo)
        end subroutine new_halo

        type(tUniverse) pure function getUniverse_halo(halo)
                type(tHalo), intent(in) :: halo
                getUniverse_halo = getUniverse(halo%nufit)
        end function getUniverse_halo

        real(D) pure function getAvgDensity_halo(halo)
                use mUniverse
                type(tHalo), intent(in) :: halo
                getAvgDensity_halo = getAvgDensity( getUniverse(halo) )
        end function getAvgDensity_halo

        type(tLPS) pure function getLPS_halo(halo)
                type(tHalo), intent(in) :: halo
                getLPS_halo = getLPS(halo%nuFit)
        end function getLPS_halo

        subroutine write_halo(halo)
                use mLPS
                type(tHalo), intent(in) :: halo
                write(6, 1) '#*************************************************'
                call write( getLPS(halo) )
                if(halo%fracSub /= 0.0_D) then
                        write(6, 2) '# Substructure With  f =', halo%fracSub
                        write(6, 2) '#                   mu =', halo%mu
                else
                        write(6, 1) '# No Substructure'
                end if
                if(halo%cDistribution) then
                        write(6, 1) '# With c Distribution'
                else
                        write(6, 1) '# No c Distribution'
                end if
                write(6, 1) '#*************************************************'
1               format(a)
2               format(a, f5.2)
        end subroutine write_halo

        subroutine refitNu(halo,a)
                type(tHalo), intent(inout) :: halo
                real(D), intent(in) :: a
                call refit(halo%nuFit, a)
        end subroutine refitNu

        real(D) function getOmega_m0_halo(halo)
                use mUniverse
                type(tHalo), intent(in) :: halo
                getOmega_m0_halo = getOmega_m0( getUniverse(halo) )
        end function getOmega_m0_halo

        real(D) pure function getFracSub(halo)
                type(tHalo), intent(in) :: halo
                getFracSub = halo%fracSub
        end function getFracSub

        real(D) pure function getM_star(halo)
                type(tHalo), intent(in) :: halo
                getM_star = halo%m_star
        end function getM_star

        real(D) pure function getRun_halo(halo)
                use mLPS
                type(tHalo), intent(in) :: halo
                getRun_halo = getRun( getLPS(halo) )
        end function getRun_halo

        real(D) pure function getH_halo(halo)
                use mUniverse
                type(tHalo), intent(in) :: halo
                getH_halo = getH( getUniverse(halo) )
        end function getH_halo

        logical pure function cDistribution(halo)
                type(tHalo), intent(in) :: halo
                cDistribution = halo%cDistribution
        end function cDistribution

        pure subroutine setFracSub(halo, fracSub)
                type(tHalo), intent(inout) :: halo
                real(D), intent(in) :: fracSub
                halo%fracSub = fracSub
        end subroutine setFracSub

        pure subroutine setCDistribution(halo)
                type(tHalo), intent(inout) :: halo
                halo%cDistribution = .true.
        end subroutine setCDistribution

        real(D) pure function getSpectralIndex_halo(halo)
                use mLPS
                type(tHalo), intent(in) :: halo
                getSpectralIndex_halo = getSpectralIndex( getLPS(halo) )
        end function getSpectralIndex_halo

        real(D) pure function getMass_nu_halo(halo)
                use mLPS
                type(tHalo), intent(in) :: halo
                getMass_nu_halo = getMass_nu( getLPS(halo) )
        end function getMass_nu_halo

        subroutine normalizeBias(halo, a)
                use mIntegration
                type(tHalo), intent(inout) :: halo
                real(D), intent(in) :: a
                type(tIntegral) :: integral

                call new(integral, nuLL(halo,a), nuUL(halo,a), &
                         eps = EPS, changeVariables = 'EXP', &
                         params = (/a/), halo = halo)
                call integrate(integral, biasNormIntegrand, halo%biasNorm)
                halo%biasNorm = 1.0_D / halo%biasNorm
                call setBiasNorm(halo, halo%biasNorm)
                call setNormalizedForA(halo, a)
        end subroutine normalizeBias

        real(D) pure function biasNorm(halo)
                type(tHalo), intent(in) :: halo
                biasNorm = halo%biasNorm
        end function biasNorm

        pure subroutine setBiasNorm(halo, biasNorm)
                type(tHalo), intent(inout) :: halo
                real(D), intent(in) :: biasNorm
                halo%biasNorm = biasNorm
        end subroutine setBiasNorm

        real(D) pure function normalizedForA(halo)
                type(tHalo), intent(in) :: halo
                normalizedForA = halo%normalizedForA
        end function normalizedForA

        pure subroutine setNormalizedForA(halo, a)
                type(tHalo), intent(inout) :: halo
                real(D), intent(in) :: a
                halo%normalizedForA = a
        end subroutine setNormalizedForA

        real(D) function nuLL(halo,a)
                type(tHalo), intent(in) :: halo
                real(D), intent(in) :: a
                nuLL = m_nu(halo, MLL, a)
        end function nuLL

        real(D) function nuUL(halo, a)
                type(tHalo), intent(in) :: halo
                real(D), intent(in) :: a
                nuUL = m_nu(halo, MUL, a)
        end function nuUL

        real(D) function nu_m_halo(halo,nu,a)
                type(tHalo), intent(in) :: halo
                real(D), intent(in) :: nu, a
                nu_m_halo = nu_m(halo%nuFit, nu, a)
        end function nu_m_halo

        real(D) function m_nu_halo(halo,m,a)
                type(tHalo), intent(in) :: halo
                real(D), intent(in) :: m, a
                m_nu_halo = m_nu(halo%nuFit, m, a)
        end function m_nu_halo

        real(D) function dnu_halo(halo, m, a)
                type(tHalo), intent(inout) :: halo
                real(D), intent(in) :: m, a
                dnu_halo = dnu(halo%nuFit, m, a)
        end function dnu_halo

        real(D) pure function Delta_c(u,a)
                use mConstants
                type(tUniverse), intent(in) :: u
                real(D), intent(in) :: a ! scale factor a=1/1+z
                Delta_c = 3.0_D * (12.0_D * PI) ** (2.0_D / 3.0_D) / 20.0_D &
                          * Delta_cFactor(u, a)
        end function Delta_c

        !*******************************************************************
        ! HaloProfile calculates u_R(k,a), the Fourier transform of the halo
        ! profile at redshift defined by a=1/1+z.
        !*******************************************************************
        real(D) function haloProfile(halo, k, a, m, cparam)
                use mConstants
                use mUniverse
                use mSineCosineIntegrals
                type(tHalo), intent(in) :: halo
                real(D), intent(in) :: k, a, m, cparam
                real(D) :: kappa, factor, &   ! misc. factors
                           si1, ci1, si2, ci2 ! sine and cosine integral values
#ifdef M_VIR
                kappa = k * m_R( universe(halo), m ) &
                        / ( cparam * Delta(a) ** ONE_THIRD )
#else ! M_VIR
                kappa = k * m_R( getUniverse(halo), m ) &
                        / (cparam * DELTA_H_TO_ONE_THIRD)
#endif ! M_VIR
                factor = kappa * (1.0_D + cparam)
                call sici(factor, si1, ci1)
                call sici(kappa, si2, ci2)
                haloProfile = ( sin(kappa) * (si1 - si2) &
                                + cos(kappa) * (ci1 - ci2) &
                                - sin(kappa * cparam) / factor ) &
                              / dimensionless_mProfile(cparam)
        end function haloProfile

        !************************************
        ! This gives the COMOVING NFW profile
        !************************************
        real(D) function NFW(halo, r, M, a, conc)
                use mConstants
                type(tHalo), intent(in) :: halo
                real(D), intent(in) :: r, M, a, conc
                real(D) :: factor
                factor = conc * r / m_R_vir(halo, M, a)
!                factor = conc * r / m_R_180(halo, M)
                NFW = 1 + factor ! temp storage
                NFW = M * (factor / r) ** 3 &
                        / ( 4.0_D * PI * factor * NFW * NFW &
                            * dimensionless_mProfile(conc) )
        end function NFW

        !**********************************************************************
        ! Calculates c_h:=r_h/r_s as a function of the halo mass M_h using an
        ! arbitrary mass definition chosen by the parameter DELTA_H. For
        ! example, if you want to use M_180, set DELTA_H=180, and pass M_180 in
        ! as M_h. If you want to use M_vir, compile with -DM_VIR and call with
        ! M_h=M_vir.
        !
        ! The method to obtain c_h is outlined in Hu & Kravtsov (2001),
        ! astro-ph/0203169.
        !**********************************************************************
        real(D) function cParam(halo, M_h, a, flag_subhalo)
                use mConstants
                type(tHalo), intent(in) :: halo
                real(D), intent(in) :: M_h, a
                logical, intent(in), optional :: flag_subhalo
                logical :: flag_pass_subhalo
#ifndef M_VIR
!                real(D), parameter :: EPS = epsilon(1.0_R4)
!                real(D), parameter :: EPS = epsilon(1.0_R4)
                real(D), parameter :: EPS = 1e-2_D
                integer, parameter :: MAX_ITER = 100
                real(D) :: delta_v, M_v, c_v, c_v_last
                integer :: i, err
                logical :: converged
#endif
                set_flag_to_pass: if( present(flag_subhalo) ) then
                        flag_pass_subhalo = flag_subhalo
                else set_flag_to_pass
                        flag_pass_subhalo = .FALSE.
                end if set_flag_to_pass
#ifdef M_VIR
                cParam = M_vir_c(halo, M_h, a, flag_pass_subhalo)
#else ! not M_VIR
#ifdef NEW_SELJAK
                cParam = M_vir_c(halo, M_h, a, flag_pass_subhalo)
#else ! not NEW_SELJAK
                converged = .false.
                c_v_last = -1.0_D
                delta_v = Delta(halo, a)
                M_v = M_h
                iterate: do i = 1, MAX_ITER
                        cParam = M_vir_c(halo, M_v, a, flag_pass_subhalo)
                        if(abs( (cParam - c_v_last) / cParam ) <= EPS) then
                                converged = .true.
!                                write(6,*) 'Converged at',i
                                exit
                        else
                                M_v = M_h * delta_v / DELTA_H &
                                    * (cParam &
                                    * HuKravtsov( DELTA_H &
                                                  / (delta_v * cParam ** 3) &
                                      * Dimensionless_mProfile(cParam) ) &
                                      ) ** 3
                                c_v_last = cParam
                                cycle
                        end if
                end do iterate

                if(.not. converged) &
                    write(0, *) 'WARNING: cParam failed to converge for M_h=', &
                                 M_h

                ! Scale cParam from r_vir/r_s to r_h/r_s
                cParam = cParam &
                         * ( M_h * delta_v / (M_v * DELTA_H) ) ** ONE_THIRD
#endif ! not NEW_SELJAK
#endif ! not M_VIR
        end function cParam

        !**************************************************************
        ! The concentration parameter as a function of the virial mass.
        !**************************************************************
        real(D) pure function M_vir_c(halo, m, a, flag_subhalo)
                type(tHalo), intent(in) :: halo
                real(D), intent(in) :: m, a
                logical, intent(in) :: flag_subhalo

                subhalo: if(flag_subhalo) then
                        ! Bullock (1999) :
                        M_vir_c = 7.5_D * a * (m / halo%m_star) ** (-0.3_D)
                else subhalo
#ifdef MAIN_SZ
                        ! Seljak (2000):
                        M_vir_c = 10.0_D * a * (m / halo%m_star) ** (-0.2_D)
#elif NEW_SELJAK
                        ! This is c_180(M_180) from Huffenberger & Seljak (2000)
                        M_vir_c = 11.0_D * a * (m / halo%m_star) ** (-0.05_D)
#else ! not MAIN_SZ not NEW_SELJAK
                        ! These are c_vir(M_vir):
                        ! Bullock (1999):
                        M_vir_c = 9.0_D * a * (m / halo%m_star) ** (-0.13_D)
                        ! Seljak (2000):
!                        M_vir_c = 10.0_D * a * (m / halo%m_star) ** (-0.2_D)
#endif ! MAIN_SZ

                end if subhalo
        end function M_vir_c

        !**********************************************************
        ! The integral of the dimensionless NFW profile (divided by
        ! rho_s) out to a factor of x times the scale radius.
        !**********************************************************
        real(D) pure function dimensionless_mProfile(x)
                real(D), intent(in) :: x
                dimensionless_mProfile = 1.0_D + x
                dimensionless_mProfile = log(dimensionless_mProfile) &
                        - x / dimensionless_mProfile
        end function dimensionless_mProfile

        !*******************************************************************
        ! A fitting formula from Hu & Kravtsov (2001), astro-ph/0203169 to
        ! calculate the inverse of x**3*dimensionless_mProfile. See comments
        ! there.
        !*******************************************************************
        real(D) pure function HuKravtsov(m)
                real(D), intent(in) :: m
                real(D), parameter :: A_1 = 0.5116_D, &
                                      A_2 = -0.4283_D, &
                                      A_3 = -3.13e-3_D, &
                                      A_4 = -3.52e-5_D
                real(D) :: p
                p = log(m)
                p = a_2 + p * (a_3 + p * a_4)
                HuKravtsov = 1.0_D / sqrt( a_1 * m ** (2.0_D * p) + 0.5625_D ) &
                             + 2.0_D * m
        end function HuKravtsov

        !*********************************************************
        ! Overdensity of collapse (a) relative to average density.
        ! This fitting formula for Delta comes from Henry (2001):
        !*********************************************************
        real(D) function Delta(halo, a)
                use mConstants
                use mUniverse
                type(tHalo), intent(in) :: halo
                real(D), intent(in) :: a
                Delta=18.0_D * PI * PI * (1.0_D + 0.4093_D * ( a * (1.0_D &
                        / getOmega_m0(halo) - 1.0_D) ** ONE_THIRD ) &
                        ** 2.71572_D)
        end function Delta

        real(D) function m_R_vir(halo, m, a)
                use mConstants
                use mUniverse
                type(tHalo), intent(in) :: halo
                real(D), intent(in) :: m, a
                m_R_vir=( 3.0_D * m / (4.0_D * PI * Delta(halo, a) &
                        * getAvgDensity(halo)) ) ** ONE_THIRD
        end function m_R_vir

        real(D) function m_r_s(halo, m, a)
                type(tHalo), intent(in) :: halo
                real(D), intent(in) :: m, a
                m_r_s = m_R_vir(halo, m, a) / cParam(halo, m, a)
        end function m_r_s

        real(D) function R_M_180(halo, r)
                use mConstants
                use mUniverse
                type(tHalo), intent(in) :: halo
                real(D), intent(in) :: r
                R_M_180 = 4.0_D / 3.0_D * PI * r ** 3 * DELTA_H &
                          * getAvgDensity(halo)
        end function R_M_180

        real(D) pure function m_r_180(halo, m)
                use mConstants
                use mUniverse
                type(tHalo), intent(in) :: halo
                real(D), intent(in) :: m
                m_r_180 = ( 3.0_D * m / (4.0_D * PI * DELTA_H * getAvgDensity(halo)) ) ** ONE_THIRD
        end function m_r_180

        !***************************************
        ! The mass function of Sheth and Tormen.
        !***************************************
        real(D) pure function STf_dnu(nu)
                real(D), intent(in) :: nu
                STf_dnu = SHETH_TORMEN_A * nu ! temp storage
                STf_dnu = SHETH_TORMEN_NORM &
                          * ( 1.0_D + (STf_dnu) ** (-SHETH_TORMEN_P) ) &
                          * sqrt(SHETH_TORMEN_A / nu) * exp(-0.5_D * STf_dnu)
        end function STf_dnu

        !**************************************************************
        ! The halo bias parameter appropriate for the Sheth-Tormen mass
        ! function.
        !**************************************************************
        real(D) pure function ST_b1(halo, nu, a)
                type(tHalo), intent(in) :: halo
                real(D), intent(in) :: nu, a
                ST_b1 = SHETH_TORMEN_A * nu ! temp storage
                ST_b1 = 1.0_D + ( ST_b1 - 1.0_D &
                      + 2.0_D * SHETH_TORMEN_P &
                              / (1.0_D + ST_b1 ** SHETH_TORMEN_P) ) &
                        / Delta_c( getUniverse(halo), a )
        end function ST_b1

        !***********************************
        ! The Press-Schechter mass function.
        !***********************************
        real(D) pure function PSf_dnu(nu)
                use mConstants
                use mLPS
                real(D), intent(in) :: nu
                PSf_dnu = sqrt( 1.0_D / (2.0_D * PI * nu) ) * exp(-0.5_D * nu)
        end function PSf_dnu

        !*****************************************************************
        ! The halo bias parameter appropriate for the Press-Schechter mass
        ! function.
        !*****************************************************************
        real(D) pure function PS_b1(halo, nu, a)
                type(tHalo), intent(in) :: halo
                real(D), intent(in) :: nu, a
                PS_b1 = 1.0_D + (nu - 1.0_D) / Delta_c( getUniverse(halo), a )
        end function PS_b1

        subroutine checkBiasNorm(halo, a)
                use mIntegration
                type(tHalo), intent(inout) :: halo
                real(D), intent(in) :: a
                type(tIntegral) :: integral
                real(D) :: norm
                if( normalizedForA(halo) /= a ) then
                        call new(integral, nuLL(halo,a), nuUL(halo,a), &
                            eps = EPS, changeVariables = 'EXP', &
                            params = (/a/), halo = halo)
                        call integrate(integral, biasNormIntegrand, norm)
                        norm = 1.0_D / norm
                        call setBiasNorm(halo, norm)
                        call setNormalizedForA(halo, a)
                end if
        end subroutine checkBiasNorm

        real(D) pure function biasNormIntegrand(integral, nu)
                use mIntegration
                type(tIntegral), intent(in) :: integral
                real(D), intent(in) :: nu
                real(D) :: a
                type(tHalo) :: halo
                a = getParam(integral)
                halo = getHalo(integral)
#ifdef PRESS_SCHECHTER
                biasNormIntegrand = PSf_Dnu(nu) * PS_b1(halo, nu, a)
#else ! PRESS_SCHECHTER
                biasNormIntegrand = STf_Dnu(nu) * ST_b1(halo, nu, a)
#endif ! PRESS_SCHECHTER
        end function biasNormIntegrand

        !*********************************************************************
        ! I divide the sub-halo mass function like this:
        ! dn/dm = M*N_0(M)*Hm(m)/m
        !
        ! N_0(M), the normalization of the sub-halo mass function. This can
        !         depend on the parent halo mass.
        !
        ! Hm(m), the sub-halo mass dependent part, with a factor of 1/m pulled
        !        out.
        !
        ! The reason for this division is that I don't want to include
        ! N_0(M) in the dm (sub-halo) integral, for performance reasons. The
        ! extra factor of M is usually cancelled by dM N(M) -> d(nu)/M f(nu).
        ! The factor of 1/m is usually cancelled by a factor of (m/rho).
        !*********************************************************************
        real(D) pure function N_0(halo,doUL)
                type(tHalo), intent(in) :: halo
                real(D), intent(in) :: doUL
#ifdef __sun
                interface erf ! interface to erf in Sun Fortran library.
                        real(4) function r_erf(x)
                                real(4), intent(in) :: x
                        end function r_erf
                        real(8) function d_erf(x)
                                real(8), intent(in) :: x
                        end function d_erf
                end interface ! erf
#endif
                real(D) :: one_minus_mu
                one_minus_mu = 1.0_D - halo%mu
                N_0 = halo%fracSub * one_minus_mu &
                        / ( doUL ** one_minus_mu &
                            - halo%subLL ** one_minus_mu )
        end function N_0

        real(D) pure function subMassFunctionHm(halo, m)
                type(tHalo), intent(in) :: halo
                real(D), intent(in) :: m
                subMassFunctionHm = m ** (-halo%mu)
        end function subMassFunctionHm

        pure subroutine SetupSubhaloLimits(halo)
                type(tHalo), intent(inout) :: halo
#ifdef MAIN_SZ
!                subLL = halo%frac_sub * 1e8_D
!                subLL = halo%frac_sub * M_min * h
!                subLL = M_min * h
!                subLL = halo%frac_sub * LL
                halo%subLL = 1e11_D
                halo%subUL = halo%frac_sub * M_max * h
!                subLL = 0.5_D * subUL
#else
!                halo%subLL = halo%frac_sub * LL
!                halo%subUL = halo%frac_sub * UL
                halo%subLL = SUBHALO_UL_FRACTION * MLL
                halo%subUL = SUBHALO_UL_FRACTION * MUL
#endif ! MAIN_SZ
        end subroutine SetupSubhaloLimits

        real(D) pure function subLL(halo)
                type(tHalo), intent(in) :: halo
                subLL = halo%subLL
        end function subLL

        real(D) pure function subUL(halo)
                type(tHalo), intent(in) :: halo
                subUL = halo%subUL
        end function subUL

        real(D) function N_gal(m)
                real(D), intent(in) :: m
                real(D), parameter :: M_blue = 4e12_D, & ! M_sun / h
                                      M_red = 2.5e12_D   ! M_sun / h
                write(0, *) 'ERROR: Factors of h are screwed up here!'
                stop
                mass: if(m < 1e11_D) then
                        N_gal = 0.0_D
                else if(m <= M_blue) then mass
                        N_gal = 0.7_D + (m / M_red) ** 0.9_D
                else mass
                        N_gal = 0.7_D * (m / M_blue) ** 0.8_D &
                                + (m / M_red) ** 0.9_D
                end if mass
        end function N_gal

        real(D) function N_sub(halo, M)
                type(tHalo), intent(in) :: halo
                real(D), intent(in) :: M ! total halo dark-matter mass
                real(D) :: minusMu, oneMinusMu, ul_sub
                minusMu = -halo%mu
                oneMinusMu = 1.0_D + minusMu
                ul_sub = SUBHALO_UL_FRACTION * M
                ! Note M_UL_FOR_G is set to subUL
                N_sub = oneMinusMu * getFracSub(halo) * M &
                        * (M_LL_FOR_G ** minusMu - M_UL_FOR_G ** minusMu) &
                        / (halo%mu &
                        * ( ul_sub ** oneMinusMu - subLL(halo) ** oneMinusMu ))
        end function N_sub

       !************************************************************
       ! The parameter describing deviation from Poisson statistics,
       ! where alpha = 1.
       !************************************************************
        real(D) function alpha(halo, M)
                use mUniverse
                type(tHalo), intent(in) :: halo
                real(D), intent(in) :: M
                if(M < 1e13_D / getH( getUniverse(halo) )) then
!                        alpha = 0.5_D * log(M * h( universe(halo) ) / 1e11_D)
                        alpha = 0.0_D
                else
                        alpha = 1.0_D
                end if
        end function alpha

        !******************************************************************
        ! Calculates the (dimensionful) power spectrum using the halo model
        !******************************************************************
        subroutine powerSpectrum(halo, k, a, PS)
                use mObject
                type(tHalo), intent(in) :: halo
                real(D), intent(in) :: k, a
                real(D), intent(out) :: PS
                real(D) :: contrib_1h, contrib_2h, &
                           contrib_1c, contrib_sc, contrib_2c

                call P_1h(halo, k, a, contrib_1h)
                call P_2h(halo, k, a, contrib_2h)
                PS = contrib_1h + contrib_2h
                substructure: if( getFracSub(halo) /= 0.0_D ) then
                        call P_1c(halo, k, a, contrib_1c)
                        call P_sc(halo, k, a, contrib_sc)
                        call P_2c(halo, k, a, contrib_2c)
                        PS = PS + contrib_1c + contrib_sc + contrib_2c
                end if substructure
        end subroutine powerSpectrum

        !****************************************************************
        ! This function calculates the one-halo term of the dimensionless 
        ! dark matter power spectrum, including substructure effects.
        ! It uses the Sheth-Tormen mass function,
        ! but this is easily changed in the function P_1hIntegrand.
        !****************************************************************
!        real(D) function P_1h(halo, k, a)
        subroutine P_1h(halo, k, a, ans)
                use mObject
                use mIntegration
                type(tHalo), intent(in) :: halo
                real(D), intent(in) :: k, a
                type(tIntegral) :: P_1hIntegral
                real(D), intent(out) :: ans
                real(D) :: one_minus_frac_sub
                call new(P_1hIntegral, nuLL(halo, a), nuUL(halo, a), &
                    eps = EPS, changeVariables = 'EXP', &
                    params = (/ k, a /), halo = halo)
                call integrate(P_1hIntegral, P_1hIntegrand, ans)

!                if( .not. converged(P_1hIntegral) ) &
!                    write(0, *) 'WARNING: P_1h failed to converge.'
                call checkConvergence(P_1hIntegral, 'P_1h failed to converge.')


                one_minus_frac_sub = 1.0_D - getFracSub(halo)
                ans = one_minus_frac_sub * one_minus_frac_sub * ans &
                    / getAvgDensity(halo)

        end subroutine P_1h

        real(D) function P_1hIntegrand(integral, nu)
                use mConstants
                use mIntegration
                use mObject
                type(tIntegral), intent(in) :: integral
                real(D), intent(in) :: nu
                real(D) :: k, a, lnCBar, u, M, M_s
                type(tHalo) :: halo
                type(tIntegral) :: P_1hCIntegral
                call getParams(integral, k, a)
                halo = getHalo(integral)

                M = nu_m(halo, nu, a)
                cDist: if( cDistribution(halo) ) then
                        lnCBar = log( CParam(halo, &
                            M * ( 1.0_D - getFracSub(halo) ), a) )
                        call new(P_1hCIntegral, &
                                 lnCBar - 3.0_D * SIGMA_LN_C, &
                                 lnCBar + 3.0_D * SIGMA_LN_C, &
                                 eps = EPS, &
                                 params = (/ k, a, M /), halo = halo )
                        call integrate(P_1hCIntegral, P_1hCIntegrand, u)
                        P_1hIntegrand = M * u &
                            / sqrt(2.0_D * PI * SIGMA_LN_C2) * STf_dnu(nu)
                        ! N.B.: there ARE two factors of u in P_1hIntegrand,
                        ! since P_1h_cIntegrand integrates u**2
                else cDist
                        M_s = M * ( 1.0_D - getFracSub(halo) )
                        u = haloProfile( halo, k, a, M_s, CParam(halo, M_s, a) )
                        P_1hIntegrand = M * u * u &
#ifdef PRESS_SCHECHTER
                            *PSf_dnu(nu)
#else ! PRESS_SCHECHTER
                            *STf_dnu(nu)
#endif ! PRESS_SCHECHTER
                end if cDist
        end function P_1hIntegrand

        real(D) function P_1hCIntegrand(integral, lnC)
                use mIntegration
                use mObject
                type(tIntegral), intent(in) :: integral
                real(D), intent(in) :: lnC
                real(D) :: k, a, M, lnCBar, u, M_s
                type(tHalo) :: halo
                
                call getParams(integral, k, a, M)
                halo = getHalo(integral)
                M_s = M * ( 1.0_D - getFracSub(halo) )
                lnCBar = log( CParam(halo, M_s, a) )
                u = haloProfile( halo, k, a, M_s, exp(lnC) )
                P_1hCIntegrand = lnC - lnCBar
                P_1hCIntegrand = u * u &
                    * exp( -P_1hCIntegrand * P_1hCIntegrand &
                           / (2.0_D * SIGMA_LN_C2) )
        end function P_1hCIntegrand

        !****************************************************************
        ! This function calculates the two-halo term of the dimensionless 
        ! dark matter power spectrum. It uses the Press-Schechter mass
        ! function and bias parameter, but this is easily changed in the
        ! function P_2hIntegrand
        !****************************************************************
!        real(D) function P_2h(halo, k, a)
        subroutine P_2h(halo, k, a, ans) 
                use mObject
                use mUniverse
                use mLPS
                use mIntegration
#ifdef NON_LINEAR_CORRELATION
                use PSFittingFunctions
#endif ! NON_LINEAR_CORRELATION
!DEREK                type(tHalo), intent(inout) :: halo
                type(tHalo), intent(in) :: halo
                real(D), intent(in) :: k, a
                type(tIntegral) :: P_2hIntegral
                real(D), intent(out) :: ans
                real(D) :: factor
!                call checkBiasNorm(halo, a)
#ifdef HALO_EXCLUSION
DEREK---broken!
                halo%contrib(2) = m_nu(LL, halo%a)
                call integrate( P_2hIntegrand, halo%contrib(2), &
                    m_nu( min(UL, R_M_180( 1.0_D / (2.0_D * k) )), halo%a ), &
                    info, .false., 'EXP', EPS, &
                    (/ k, halo%a, halo%cDistribution /) )
!#elif NON_LINEAR_CORRELATION
!DEREK HERE
#else ! not HALO_EXCLUSION, not NON_LINEAR_CORRELATION

!                if( cDistribution(halo) ) &
!                    write(0, *) 'P_2h w/ c-dist not right!'

                call new(P_2hIntegral, nuLL(halo, a), nuUL(halo, a), &
                         eps = EPS, changeVariables = 'EXP', &
                         params = (/ k, a /), halo = halo )
                call integrate(P_2hIntegral, P_2hIntegrand, ans)
#endif ! not HALO_EXCLUSION, not NON_LINEAR_CORRELATION

!                if( .not. converged(P_2hIntegral) ) &
!                    write(0, *) 'WARNING: P_2h failed to converge.'
                call checkConvergence(P_2hIntegral, 'P_2h failed to converge.')

                ans = ans * biasNorm(halo)
#ifdef NON_LINEAR_CORRELATION
DEREK---broken
                halo%contrib(2)=2.0_d*PI*PI/k**3*SmithEtAl(k,halo%a) &
                    *halo%contrib(2)*halo%contrib(2)
#else ! not NON_LINEAR_CORRELATION
                factor = growthFactor( getUniverse(halo), a ) * ans
                ans = factor * factor * P_L( getLPS(halo), k, a )
#endif ! NON_LINEAR_CORRELATION
        end subroutine P_2h

        real(D) function P_2hIntegrand(integral, nu)

                use mConstants
                use mIntegration
                use mObject

                type(tIntegral), intent(in) :: integral
                real(D), intent(in) :: nu
                real(D) :: k, a, parentM, lnCBar, u
                type(tIntegral) :: P_2h_cIntegral
                type(tHalo) :: halo

                call getParams(integral, k, a)
                halo = getHalo(integral)

                parentM = nu_m(halo, nu, a)

                if(cDistribution( getHalo(integral)) ) then
                        lnCBar = log( &
                            CParam( getHalo(integral), parentM, a, .false. ))
                        call new(P_2h_cIntegral, &
                                 lnCBar - 3.0_D * SIGMA_LN_C, &
                                 lnCBar + 3.0_D * SIGMA_LN_C, &
                                 eps = EPS, &
                                 params = (/ k, a, parentM /), &
                                 halo = getHalo(integral))
                        call Integrate(P_2h_cIntegral, P_2hCIntegrand, &
                                       P_2hIntegrand)

                        P_2hIntegrand = P_2hIntegrand &
                                        * ST_b1( getHalo(integral), nu, a ) &
                                        / sqrt( 2.0_D * PI * SIGMA_LN_C2) &
                                        * STf_dnu(nu)
                else
#ifdef NON_LINEAR_CORRELATION
DEREK---broken!
                P_2hIntegrand=(1.0_d+0.2_d &
                    *Smith2Point(1.0_d/kac(1),kac(2)))**(-0.25_d) &
                    *ST_b1(nu,kac(2)) &
                    *haloProfile(kac(1),kac(2), parentM, 
                    ,CParam(P_2hIntegrand,kac(2))) &
                    *STf_dnu(nu)
#else ! NON_LINEAR_CORRELATION
#ifdef PRESS_SCHECHTER
                P_2hIntegrand = PSf_dnu(nu) * PS_b1(halo, nu, a) &
#else ! PRESS_SCHECHTER
                P_2hIntegrand = STf_dnu(nu) * ST_b1(halo, nu, a) &
#endif ! PRESS_SCHECHTER
                    * haloProfile( halo, k, a, parentM, &
                        cParam(halo, parentM, a) )
#endif ! NON_LINEAR_CORRELATION
                end if
        end function P_2hIntegrand

        real(D) function P_2hCIntegrand(integral, lnC)
                use mIntegration
                use mObject
                type(tIntegral), intent(in) :: integral
                real(D), intent(in) :: lnC
!                real(D), intent(in), optional :: kaM(:) ! kaM(1) = k
!                                                        ! kaM(2) = a
!                                                        ! kaM(3) = M
                real(D) :: k, a, M, lnCBar, u

                call getParams(integral, k, a, M)

                lnCBar = cParam( getHalo(integral), M, a, .false. )
                u = haloProfile( getHalo(integral), k, a, M, exp(lnC) )
                lnCBar = log(lnCBar)
                P_2hCIntegrand = lnC - lnCBar
                P_2hCIntegrand = u * exp( -P_2hCIntegrand * P_2hCIntegrand &
                                           / (2.0_D * SIGMA_LN_C2) )
        end function P_2hCIntegrand

#ifdef NON_LINEAR_CORRELATION
DEREK---FFT is broken
        real(D) function smith2Point(r,a)
                use mFFT
                real(D), intent(in) :: r, a
                real(D), parameter :: FTkLL = 1e-2_D, &
                                      FTkUL=1000.0_D
                REAl(D) :: trash ! Not used
                call FFTIntegrate( smith2PointIntegrand, FTkLL, FTkUL, &
                                   r, Smith2Point, trash, (/ a /) )
                smith2Point = smith2Point / r
        end function smith2Point

        real(D) function smith2PointIntegrand(k, a)
                use mPSFittingFunctions
                real(D), intent(in) :: k
                real(D), intent(in), optional :: a(:)
                smith2PointIntegrand = smithEtAl( k, a(1) ) / (k * k)
        end function smith2PointIntegrand
#endif ! NON_LINEAR_CORRELATION

        ! See Sheth & Jain, astro-ph/0208353
        subroutine P_1c(halo, k, a, ans)
                use mObject
                use mIntegration
                type(tHalo), intent(in) :: halo
                real(D), intent(in) :: k, a
                type(tIntegral) :: P_1cIntegral
                real(D), intent(out) :: ans

                ! I define the lower limit of parent halos as
                ! subLL / SUBHALO_UL_FRACTION because parent halos
                ! smaller than this should not contain sub-halos.
                call new(P_1cIntegral, &
                         m_nu( halo, subLL(halo) / SUBHALO_UL_FRACTION, a ), &
                         nuUL(halo, a), eps = EPS, changeVariables = 'EXP', &
                         params = (/ k, a /), halo = halo)
                call integrate(P_1cIntegral, P_1cIntegrand, ans)

!                if( .not. converged(P_1cIntegral) ) then
!                        write(0, *) 'WARNING: P_1cIntegrand integration failed'
!                        write(0, *) 'to converge.'
!                end if
                call checkConvergence(P_1cIntegral, 'P_1cIntegrand integration failed to converge.')

                ans = ans / getAvgDensity(halo)
        end subroutine P_1c

        real(D) function P_1cIntegrand(integrand, nu)
                use mIntegration
                use mObject
                type(tIntegral), intent(in) :: integrand
                real(D), intent(in) :: nu
                type(tIntegral) :: P_1cSubIntegral
                real(D) :: k, a, parentM, doUL
                type(tHalo) :: halo

                call getParams(integrand, k, a)
                halo = getHalo(integrand)

                parentM = nu_m(halo, nu, a)
                doUL = min( subUL(halo), SUBHALO_UL_FRACTION * parentM)
                call new(P_1cSubIntegral, subLL(halo), doUL, eps = EPS, &
                         changeVariables = 'EXP', &
                         params = (/ k, a /), halo = halo )
                call integrate(P_1cSubIntegral, P_1cSubIntegrand, P_1cIntegrand)

!                if( .not. converged(P_1cSubIntegral) ) &
!                    write(0, *) 'WARNING: P_cSubIntegrand failed to converge.'
                call checkConvergence(P_1cSubIntegral, 'P_cSubIntegrand failed to converge.')

                P_1cIntegrand = N_0(halo, doUL) * P_1cIntegrand
                P_1cIntegrand = P_1cIntegrand * STf_dnu(nu)
        end function P_1cIntegrand

        real(D) function P_1cSubIntegrand(integral, m)
                use mIntegration
                use mObject
                type(tIntegral), intent(in) :: integral
                real(D), intent(in) :: m
                real(D) :: k, a, u
                type(tHalo) :: halo

                call getParams(integral, k, a)
                halo = getHalo(integral)

                u = haloProfile( halo, k, a, m, cParam(halo, m, a, .true.) )
                P_1cSubIntegrand = subMassFunctionHm(halo, m) * m * u * u
        end function P_1cSubIntegrand

        subroutine P_sc(halo, k, a, ans)
                use mIntegration
                use mObject
                type(tHalo), intent(in) :: halo
                real(D), intent(in) :: k, a
                type(tIntegral) :: P_scIntegral
                real(D), intent(out) :: ans

                ! I define the lower limit of parent halos as
                ! subLL / SUBHALO_UL_FRACTION because halos smaller than that
                ! should not contain sub-halos.
                call new(P_scIntegral, &
                         m_nu( halo, subLL(halo) / SUBHALO_UL_FRACTION, a ), &
                         nuUL(halo, a), eps = EPS, changeVariables = 'EXP', &
                         params = (/ k, a /), halo = halo)
                call integrate(P_scIntegral, P_scIntegrand, ans)

!                if( .not. converged(P_scIntegral) ) &
!                    write(0, *) 'WARNING: P_scIntegrand failed to converge.'
                call checkConvergence(P_scIntegral, 'P_scIntegrand failed to converge.')

                ans = 2.0_D * ( 1.0_D - getFracSub(halo) ) &
                    * ans / getAvgDensity(halo)
        end subroutine P_sc

        real(D) function P_scIntegrand(integral, nu)
                use mIntegration
                use mObject
                type(tIntegral), intent(in) :: integral
                real(D), intent(in) :: nu
                real(D) :: k, a
                type(tHalo) :: halo
                real(D) :: parentU, parentM, doUL, M_s
                type(tIntegral) :: P_scSubIntegral

                call getParams(integral, k, a)
                halo = getHalo(integral)

                parentM = nu_m(halo, nu, a)
                doUL = min( subUL(halo), SUBHALO_UL_FRACTION * parentM)
                call new(P_scSubIntegral, subLL(halo), doUL, eps = EPS, &
                    changeVariables = 'EXP', &
                    params = (/ k, a /), halo = halo)
                call integrate(P_scSubIntegral, P_scSubIntegrand, P_scIntegrand)

!                if( .not. converged(P_scSubIntegral) ) &
!                    write(0, *) 'WARNING: P_scSubIntegrand failed to converge.'
                call checkConvergence(P_scSubIntegral, 'P_scSubIntegrand failed to converge.')

                P_scIntegrand = N_0(halo, doUL) * parentM * P_scIntegrand
                M_s = parentM * ( 1.0_D - getFracSub(halo) )
                parentU = haloProfile( halo, k, a, M_s, &
                    cParam(halo, M_s, a, .false.) )
                P_scIntegrand = STf_dnu(nu) * parentU * parentU * P_scIntegrand
        end function P_scIntegrand

        real(D) function P_scSubIntegrand(integral, m)
                use mIntegration
                use mObject
                type(tIntegral), intent(in) :: integral
                real(D), intent(in) :: m
                real(D) :: k, a
                type(tHalo) :: halo

                call getParams(integral, k, a)
                halo = getHalo(integral)

                P_scSubIntegrand = subMassFunctionHm(halo, m) &
                    * haloProfile( halo, k, a, m, cParam(halo, m, a, .true.) )

        end function P_scSubIntegrand

!        real(D) function P_2c(halo, k, a)
        subroutine P_2c(halo, k, a, ans)
                use mIntegration
                use mObject
                type(tHalo), intent(in) :: halo
                real(D), intent(in) :: k, a
                type(tIntegral) :: P_2cIntegral
                real(D), intent(out) :: ans

                ! I define the lower limit of parent halos as
                ! subLL / SUBHALO_UL_FRACTION because halos smaller than that
                ! should not contain sub-halos.
                call new(P_2cIntegral, &
                         m_nu( halo, subLL(halo) / SUBHALO_UL_FRACTION, a), &
                         nuUL(halo,a), &
                         eps = EPS, changeVariables='EXP', &
                         params = (/ k, a /), halo = halo)
                call integrate(P_2cIntegral, P_2cIntegrand, ans)

!                if( .not. converged(P_2cIntegral) ) &
!                    write(0, *) 'WARNING: P_2cIntegrand failed to converge.'
                call checkConvergence(P_2cIntegral, 'P_2cIntegrand failed to converge.')

                ans = ans / getAvgDensity(halo)
        end subroutine P_2c

        real(D) function P_2cIntegrand(integral, nu)
                use mIntegration
                use mObject
                type(tIntegral), intent(in) :: integral
                real(D), intent(in) :: nu
                real(D) :: k, a, doUL, parentM, M_s
                type(tHalo) :: halo
                type(tIntegral) :: P_2cSubIntegral

                call getParams(integral, k, a)
                halo = getHalo(integral)

                parentM = nu_m(halo, nu, a)
                doUL = min( subUL(halo), SUBHALO_UL_FRACTION * parentM)
                call new(P_2cSubIntegral, subLL(halo), doUL, &
                         eps = EPS, changeVariables = 'EXP', &
                         params = (/ k, a /), halo = halo )
                call integrate(P_2cSubIntegral, P_2cSubIntegrand, P_2cIntegrand)

!                if(.not.converged(P_2cSubIntegral)) &
!                    write(0,*)'WARNING: P_ccSubIntegrand failed to converge.'
                call checkConvergence(P_2cSubIntegral, 'P_ccSubIntegrand failed to converge.')

                P_2cIntegrand = N_0(halo, doUL) * P_2cIntegrand
                M_s = parentM * ( 1.0_D - getFracSub(halo) )
                P_2cIntegrand = P_2cIntegrand * haloProfile( halo, k, a, M_s, &
                    CParam(halo, M_s, a, .false.) )
                P_2cIntegrand = STf_dnu(nu) * parentM &
                                * P_2cIntegrand * P_2cIntegrand
        end function P_2cIntegrand

        real(D) function P_2cSubIntegrand(integral, m)
                use mIntegration
                use mObject
                type(tIntegral), intent(in) :: integral
                real(D), intent(in) :: m
                real(D) :: k, a
                type(tHalo) :: halo

                call getParams(integral, k, a)
                halo = getHalo(integral)

                P_2cSubIntegrand = subMassFunctionHm(halo, m) &
                    * haloProfile( halo, k, a, m, CParam(halo, m, a, .true.) )
        end function P_2cSubIntegrand

end module mHaloModel
